CITY_US
I Предварительный анализ данных
if (interactive() && !str_ends(getwd(), "R/Statistical_Analysis_2022/CITY_US/Tkachenko")) {
setwd("R/Statistical_Analysis_2022/CITY_US/Tkachenko")
}
data <- read_excel("../CITY_US.STD/CITY_shortname.xls")
data[data == "NA"] <- NA
data[, -(1:2)] <- data.frame(lapply(data[, -(1:2)], as.numeric))
fullnames <- names(read_excel("../CITY_US.STD/CITY.xls"))1 Разобраться в том, что означают признаки.
print(fullnames)## [1] "CITY City"
## [2] "STATE State"
## [3] "AREA Land area, 1990 (sq.miles)"
## [4] "R_AREA Land area, 1990 (sq.miles), ranks"
## [5] "POP92 Population,1992"
## [6] "R_POP92 Population,1992, ranks"
## [7] "POP80 Population,1980"
## [8] "R_POP80 Population,1980, ranks"
## [9] "POP_CH Population growth rate,1980-1992"
## [10] "R_POP_CH Population growth rate,1980-1992, ranks"
## [11] "POPDEN Population per square mile, 1992"
## [12] "R_POPDEN Population per square mile, 1992, ranks"
## [13] "OLD % older 65 years"
## [14] "R_OLD % older 65 years, ranks"
## [15] "BLACK Black population,1990"
## [16] "R_BLACK Black population,1990, ranks"
## [17] "BLACK% % Black population,1990"
## [18] "R_BLACK% % Black population,1990, ranks"
## [19] "ASIAN Asian or Pacific Islander population,1990"
## [20] "R_ASIAN Asian or Pacific Islander population,1990, ranks"
## [21] "ASIAN% % Asian or Pacific Islander population,1990"
## [22] "R_ASIAN% % Asian or Pacific Islander population,1990, ranks"
## [23] "HISP Hispanic population, 1990"
## [24] "R_HISP Hispanic population, 1990, ranks"
## [25] "HISP% % Hispanic population, 1990"
## [26] "R_HISP% % Hispanic population, 1990, ranks"
## [27] "BORN_F % persons of foreign born"
## [28] "R_BORN_F % persons of foreign born, ranks"
## [29] "LANG % of persons, speaking language other than English at home, 1990"
## [30] "R_LANG % of persons, speaking language other than English at home, 1990, ranks"
## [31] "HH1 % 1-person households, 1990"
## [32] "R_HH1 % 1-person households, 1990, ranks"
## [33] "FAMIL1 % one-parent family households, 1990"
## [34] "R_FAMIL1 % one-parent family households, 1990, ranks"
## [35] "DEATH Infant death rate per 1000 live births,1988"
## [36] "R_DEATH Infant death rate per 1000 live births,1988, ranks"
## [37] "CRIME Crime rate per 100000 population, 1991"
## [38] "R_CRIME Crime rate per 100000 population, 1991, rate"
## [39] "SCHOOL % of elementary and high school enrollment in public schools, 1990"
## [40] "R_SCHOOL % of elementary and high school enrollment in public schools, 1990, ranks"
## [41] "DEGREE % of adults with a bachelor's degree or higher, 1990"
## [42] "R_DEGREE % of adults with a bachelor's degree or higher, 1990, ranks"
## [43] "INCOME Median household income, 1989"
## [44] "R_INCOME Median household income, 1989, ranks"
## [45] "ASSIST % of households, receiving public assistance,1989"
## [46] "R_ASSIST % of households, receiving public assistance,1989, ranks"
## [47] "POVERT % of persons below poverty level, 1989"
## [48] "R_POVERT % of persons below poverty level, 1989, ranks"
## [49] "OLB_BIL % of housing units built 1939 or earlier, 1990"
## [50] "R_OLD_BI % of housing units built 1939 or earlier, 1990, ranks"
## [51] "OWNER Median value of specified owner-occupied housing units ($), 1990"
## [52] "R_OWNER Median value of specified owner-occupied housing units ($), 1990, ranks"
## [53] "RENTER % renter-occupied housing units,1990"
## [54] "R_RENTER % renter-occupied housing units,1990, ranks"
## [55] "GROSS Median gross rent of specified renter-occupied housing units ($), 1990"
## [56] "R_GROSS Median gross rent of specified renter-occupied housing units ($), 1990, ranks"
## [57] "CONDOM % condominimum occupied housing units, 1990"
## [58] "R_CONDOM % condominimum occupied housing units, 1990, ranks"
## [59] "TRANSP % of workers using public transportation"
## [60] "R_TRANSP % of workers using public transportation, ranks"
## [61] "UNEMP Unemployment rate, 1991"
## [62] "R_UNEMP Unemployment rate, 1991, ranks"
## [63] "LABOR Labor force - % change 1980-1990"
## [64] "R_LABOR Labor force - % change 1980-1990, ranks"
## [65] "LAB_F Female civilian labor force participation rate, 1990"
## [66] "R_LAB_F Female civilian labor force participation rate, 1990, ranks"
## [67] "MANLAB % of civilian labor force emploied in manufacturing, 1990"
## [68] "R_MANLAB % of civilian labor force emploied in manufacturing, 1990, ranks"
## [69] "TAXE City government taxes per capita, 1991"
## [70] "R_TAXE City government taxes per capita, 1991, ranks"
## [71] "TEMPER Average daily temperature in July"
## [72] "R_TEMPER Average daily temperature in July, ranks"
## [73] "PRECEP Annual precipitation"
## [74] "R_PRECEP Annual precipitation"2 Отобрать признаки
names_interesting <- c("AREA", "POP80", "POP92", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")
data <- data %>% select(all_of(c("CITY", "STATE", names_interesting)))
print(head(data))## # A tibble: 6 × 12
## CITY STATE AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 NEW_… NY 309. 7.07e6 7.31e6 23671 9236 28.4 19.3 29823 8.6 76.8
## 2 LOS_… CA 469. 2.97e6 3.49e6 7436 9730 38.4 18.9 30925 9 74.3
## 3 CHIC… IL 227. 3.01e6 2.77e6 12185 NA 16.9 21.6 26301 8.4 75.1
## 4 HOUS… TX 540. 1.60e6 1.69e6 3131 10824 17.8 20.7 26261 6.1 83.5
## 5 PHIL… PA 135. 1.69e6 1.55e6 11492 6835 6.6 20.3 24603 8 76.7
## 6 SAN_… CA 324 8.76e5 1.15e6 3546 8537 20.9 13.4 10 6.2 713 Определить вид признаков
Город и штат качественные, остальные количественные, ранги были порядковыми.
find_mode_freq <- function(x) {
x <- x[!is.na(x)]
return(max(tabulate(match(x, x))))
}
print(data %>% summarise(across(
all_of(names_interesting),
find_mode_freq
)))## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <int> <int> <int> <int> <int> <int> <int> <int> <int> <int>
## 1 1 1 1 2 1 4 2 1 5 3print(sort(data$UNEMP))## [1] 2.3 3.2 3.8 3.9 4.1 4.4 4.6 4.7 4.7 4.8 4.8 5.0 5.0 5.0 5.0
## [16] 5.1 5.1 5.1 5.3 5.4 5.4 5.4 5.4 5.4 5.6 5.7 5.9 5.9 5.9 6.0
## [31] 6.0 6.1 6.1 6.1 6.1 6.2 6.2 6.4 6.4 6.5 6.7 6.7 6.7 6.8 6.9
## [46] 6.9 7.0 7.1 7.2 7.2 7.3 7.3 7.3 7.5 7.6 7.7 7.7 7.7 8.0 8.1
## [61] 8.4 8.4 8.5 8.6 9.0 9.0 9.4 9.5 9.6 10.0 10.4 10.6 10.7 11.0 11.9
## [76] 12.6 13.1Все количественные буду считать непрерывными. возможно UNEMP непрерывный с плохой точностью.
4 не актуально
5 Построить matrix plot
if (interactive()) pdf("ggpairs_unedited.pdf")
ggpairs(
data[, -(1:2)],
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
diag = list(continuous = "barDiag")
)
if (interactive()) dev.off()7 outliers
Убираю outliers:
Помечаю некорректные данные в INCOME как NA. Удаляю город из Аляски за плотность населения. Флорида выделсется на BORN_F-INCOME Гаваи выделяются низким уровнем безработицы. странно?
data$INCOME[data$INCOME < 100] <- NA
data <- data %>% filter(STATE != "AK")6 Несимметричные распределения
Функция, которая логарифмирует, если это сделает выборку симметричнее
log_asymmetric <- function(x) {
if (skewness(x, na.rm = TRUE) < abs(skewness(log(x), na.rm = TRUE))) {
print("default")
return(x)
} else {
print("logged")
return(log(x))
}
}Автоматически логарифмирую то что имеет асимметрию и длинный хвост справа
data_logged <- data %>%
mutate(across(all_of(names_interesting), log_asymmetric))## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "logged"
## [1] "default"
## [1] "default"
## [1] "logged"
## [1] "default"if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
data_logged[, -(1:2)],
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
diag = list(continuous = "barDiag")
)
if (interactive()) dev.off()8 однородность
выглядит однородно.
9 не актуально
10 всякие характеристики
print_characteristics <- function(x) {
list(
mean = mean(x, na.rm = TRUE),
var = var(x, na.rm = TRUE),
skewness = skewness(x, na.rm = TRUE),
kurtosis = kurtosis(x, na.rm = TRUE)
)
}
sapply(data_logged %>% select(all_of
(names_interesting)), print_characteristics)## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT
## mean 4.754638 12.9069 13.01226 8.257586 9.206558 1.959747 18.11842
## var 0.6658984 0.5142076 0.4464288 0.5099096 0.06798296 0.8435754 34.68872
## skewness 0.1317646 1.453879 1.743773 0.1015321 0.1476689 0.3081197 0.2243883
## kurtosis 2.609264 6.033026 6.885752 2.846858 3.088721 2.256848 2.646043
## INCOME UNEMP TEMPER
## mean 25282.48 1.873674 77.60132
## var 14321837 0.09967048 37.59853
## skewness -0.2196584 -0.2186434 -0.1426746
## kurtosis 2.353211 3.659464 3.747206II О виде распределений и о сравнении распределений
1 выполнить первое второе задание в логичном порядке
2 анализ вида распределения признаков
2. Описать распределения признаков визуально (с помощью нормальной бумаги и PP-plot) и по критерию хи-квадрат и другим критериям.
Сюда входит: normal probability plot (что это такое?), проверка по критериям Лиллиефорса, AD, хи-квадрат, Шапиро-Уилка. По критерию хи-квадрат, а также визуально по PP-plot можно проверить и гипотезы о согласии с другими распределениями, например, логнормальным.
Рассматриваем нормальность
все плотности.
if (interactive()) pdf(file = "all_density_plots.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(value)) +
theme_bw() +
geom_histogram(aes(y = ..density..), bins = 15) +
# geom_density() +
facet_wrap(~variable, scales = "free") +
geom_line(aes(y = dnorm(value,
mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
)), color = "red") +
theme(legend.position = "bottom") +
labs(x = "", y = "")
if (interactive()) dev.off()Normal probability plot
if (interactive()) pdf(file = "normal_probability_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
stat_qq_point(size = 2) +
geom_abline() +
facet_wrap(~variable, scales = "free")
if (interactive()) dev.off()PP-plot
if (interactive()) pdf(file = "PP_plot.pdf")
ggplot(melt(data_logged, id = c("CITY", "STATE")), aes(sample = value)) +
stat_pp_point(size = 2) +
facet_wrap(~variable, scales = "free") +
stat_pp_line()
if (interactive()) dev.off()Всякие тесты на нормальность
test_all <- function(data, test, ...) {
print(test(1:10, ...)$method)
print(data %>% summarise(across(
everything(),
function(x) {
test(x, ...)$p.value
}
)))
}
test_all(data_logged[, -(1:2)], shapiro.test)## [1] "Shapiro-Wilk normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.757 0.0000113 0.000000310 0.673 0.725 0.0892 0.675 0.463 0.487 0.368test_all(data_logged[, -(1:2)], lillie.test)## [1] "Lilliefors (Kolmogorov-Smirnov) normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.536 0.0118 0.000633 0.291 0.580 0.148 0.646 0.866 0.746 0.766test_all(data_logged[, -(1:2)], ad.test)## [1] "Anderson-Darling normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.492 0.000167 0.000000752 0.370 0.620 0.0362 0.830 0.692 0.697 0.371test_all(data_logged[, -(1:2)], pearson.test)## [1] "Pearson chi-square normality test"
## # A tibble: 1 × 10
## AREA POP80 POP92 POPDEN CRIME BORN_F POVERT INCOME UNEMP TEMPER
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0.697 0.0377 0.0179 0.664 0.728 0.265 0.729 0.857 0.897 0.0626итак, все кроме населения распределено нормально (или логнормально)
1. Описать разницу между численностью населения в 1980 и 1992 годах визуально (с помощью ящиков с усами) и по критериям. Добавить признак, который разбивает города на “холодные” и “жаркие”. Сравнить их по характеристикам.
добавим признак
temper_split <- mean(data_logged$TEMPER)
data_logged$F_TEMPER <- factor(ifelse(data_logged$TEMPER < temper_split, "cold", "hot"))3 Сначала имеет смысл посмотреть на сравнение распределений в группах с помощью ящиков с усами
if (interactive()) pdf("boxplot_pop.pdf")
ggplot(
melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("POP80", "POP92")),
aes(x = variable, y = value, color = F_TEMPER)
) +
geom_boxplot() +
scale_color_manual(values = c("blue", "red"))
if (interactive()) dev.off()
if (interactive()) pdf("boxplot_pop_all.pdf")
ggplot(
melt(data_logged, id.vars = "F_TEMPER", measure.vars = c("AREA", "POP80", "POPDEN", "CRIME", "BORN_F", "POVERT", "INCOME", "UNEMP", "TEMPER")),
aes(x = variable, y = value, color = F_TEMPER)
) +
geom_boxplot() +
facet_wrap(~variable, scales = "free") +
scale_color_manual(values = c("blue", "red"))## Warning: Removed 11 rows containing non-finite values (`stat_boxplot()`).
if (interactive()) dev.off()Население растет, правый ящик с усами выше.
4 начинать нужно с t-критерия…
leveneTest(value ~ variable, data = melt(data_logged[c("POP80", "POP92")]), center = "mean")## No id variables; using all as measure variables## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 1 0.3172 0.5741
## 150t.test(
data_logged$POP80,
data_logged$POP92,
var.equal = TRUE,
paired = TRUE
)##
## Paired t-test
##
## data: data_logged$POP80 and data_logged$POP92
## t = -4.7135, df = 75, p-value = 1.098e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.14988791 -0.06083053
## sample estimates:
## mean difference
## -0.1053592Мы не доверяем критерию фишера потому что он работает только для нормальных. применим тест левена, он не отвергает гипотезу что дисперсии равны, как и ожидалось для почти одинаковых распредлений. считаем что дисперсии равны и применяем t-test.
5 если бы выборка было маленькая или дисперсии были неравны, пришлось бы применять непараметрический t-test Вилксона.
wilcox.test(data_logged$POP80, data_logged$POP92, paired = TRUE)##
## Wilcoxon signed rank test with continuity correction
##
## data: data_logged$POP80 and data_logged$POP92
## V = 730, p-value = 0.0001492
## alternative hypothesis: true location shift is not equal to 06 Смотрите на результаты применения критерия Манна-Уитни
Манна-Уитни = Вилксона, этому тесту не нужны дополонительные требования, но он менее мощен чем t-test.
7 критерий Колмогорова-Смирнова
еще есть тест колмогорова смирнова, который сравнивает распеределения.
ks.test(data_logged$POP80, data_logged$POP92)##
## Exact two-sample Kolmogorov-Smirnov test
##
## data: data_logged$POP80 and data_logged$POP92
## D = 0.15789, p-value = 0.3012
## alternative hypothesis: two-sidedон не отвергает.
8 сравнение зависимых выборок
это все и было сравнение зависимых выборок
из-за того что тесты парные, они смогли отвергнуть гипотезы о том, что средние равны, без этого гипотезы не отвеграются. статистика зависимого критерия больше статистики независимого (при положительной корреляции), поэтому мощнее.
3 Об анализе зависимостей
1 вспомним pairs plot
if (interactive()) pdf("ggpairs_logged.pdf")
ggpairs(
data_logged[, -(1:2)],
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.3)),
diag = list(continuous = "barDiag")
)
if (interactive()) dev.off()if (interactive()) pdf("ggpairs_logged_colored.pdf")
ggpairs(
data_logged[, -(1:2)],
aes(color = F_TEMPER),
lower = list(continuous = wrap("points", alpha = 0.5, size = 0.6)),
diag = list(continuous = "barDiag"),
columns = names_interesting
) +
scale_color_manual(values = c("blue", "red")) +
scale_fill_manual(values = c("blue", "red"))
if (interactive()) dev.off()Матрица корреляций пирсона
data_logged <- data_logged %>%
relocate(TEMPER, AREA, POP92, POP80, POPDEN, BORN_F, UNEMP, POVERT, CRIME, INCOME)
if (interactive()) pdf(file = "cor_matrix_pearson.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "pearson", use = "pairwise.complete.obs"))) +
geom_raster(aes(x = Var2, y = Var1, fill = value)) +
geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
scale_fill_gradient2() +
theme_dark() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
if (interactive()) dev.off()Матрица корреляций спирмана
if (interactive()) pdf(file = "cor_matrix_spearman.pdf")
ggplot(melt(cor(data_logged %>% select(-CITY, -STATE, -F_TEMPER), method = "spearman", use = "pairwise.complete.obs"))) +
geom_raster(aes(x = Var2, y = Var1, fill = value)) +
geom_text(aes(x = Var2, y = Var1, label = round(value, 2))) +
scale_fill_gradient2() +
theme_dark() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
if (interactive()) dev.off()if (interactive()) pdf(file = "all_density_groupped_plots.pdf")
ggplot(
melt(data_logged, id = c("CITY", "STATE", "F_TEMPER")),
aes(value, color = F_TEMPER)
) +
theme_bw() +
# geom_histogram(bins = 15) +
geom_density() +
facet_wrap(~variable, scales = "free") +
geom_line(aes(y = dnorm(value,
mean = tapply(value, variable, mean, na.rm = TRUE)[PANEL],
sd = tapply(value, variable, sd, na.rm = TRUE)[PANEL]
)), color = "black", linetype = 2) +
theme(legend.position = "bottom") +
labs(x = "", y = "") +
scale_color_manual(values = c("blue", "red"))
if (interactive()) dev.off()
data_logged$F_TEMPER <- NULL
# plot_all_density(data_logged_melt, aes(value, color = F_TEMPER)) +
# scale_color_manual(values = c("blue", "red"))cor(data_logged$POP80, data_logged$POP92, method = "pearson")## [1] 0.9628711cor(data_logged$POP80, data_logged$POP92, method = "spearman")## [1] 0.9155161